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ABSTRACT 

^ ■ We present a new computation of the linear tidal interaction of a protoplanetary 

Q ' core with a thin gaseous disc in which it is fully embedded. For the first time a 

discussion of the orbital evolution of cores on eccentric orbits with eccentricity (e) 
significantly larger than the gas-disc scale height to radius ratio {H/r) is given. We 
CO I find that the direction of orbital migration reverses for e > l.lH/r. This occurs as a 

. result of the orbital crossing of resonances in the disc that do not overlap the orbit 

when the eccentricity is very small. In that case resonances always give a net torque 
corresponding to inward migration. Simple expressions giving approximate fits to the 
eccentricity damping rate and the orbital migration rate are presented. We go on to 
' calculate the rate of increase of the mean eccentricity for a system of protoplanetary 

, cores due to dynamical relaxation. By equating the eccentricity damping time-scale 

with the dynamical relaxation time-scale we deduce that, for parameters thought to 
be applicable to protoplanetary discs, an equilibrium between eccentricity damping 
On ■ and excitation through scattering is attained on a 10'^-10''yr time-scale, at 1 au. This 

On ' equilibrium is maintained during the further migrational and coUisional evolution of 

fH I the system, which occurs on much longer time-scales. The equilibrium thickness of 

the protoplanet distribution is related to the equilibrium eccentricity and is such that 
it is generally well confined within the gas disc. By use of a three dimensional direct 
summation N-body code we simulate the evolution of a system of protoplanetary 
cores, initialised with a uniform isolation mass of O.IMq, incorporating our eccentricity 
damping and migration rates. Assuming that collisions lead to agglomeration, we find 
that the vertical confinement of the protoplanet distribution permits cores to build up 
in mass by a factor of ^ 10 in only ~ 10"* yr, within 1 au. The time-scale required to 
achieve this is comparable to the migration time-scale. In the context of our model and 
its particular initial conditions we deduce that it is not possible to build up a massive 
' enough core to form a gas giant planet, before orbital migration ultimately results in 

the preferential delivery of all such bodies to the neighbourhood of the central star. 
It remains to be investigated whether different disc models or initial planetesimal 
distributions might be more favourable for slowing or halting the migration, leading 
to possible giant planet formation at intermediate radii. 

Key words: accretion, accretion discs - celestial mechanics, stellar dynamics - Solar 
system: formation - stars: formation - planetary systems 
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1 INTRODUCTION 

The last decade has seen the discovery of Jupiter-sized extrasolar planets orbiting Solar-type stars. In roughly one fifth of the 
cases recorded to date a planet is inferred to be in nearly circular orbit (having eccentricities less than 0.05) at a distance 
from the parent star of less than fifteen Solar radii (~ 0.06 au). In the rest of the detections the planets are inferred to have 
eccentricities up to 0.7 (with a mean value of 0.3), for semimajor axes in the range 0.07-3 au (Vogt et al. 2000). Thus the 
orbital characteristics of many of the known extrasolar systems contrast strongly with the gas giants of the Solar system 
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that have eccentricities less than about 0.05 and semimajor axes greater than 5au. These new discoveries pose a considerable 
challenge to the standard models of planetary formation that are based on our knowledge of the Solar system alone. 

According to the accumulation model, protoplanetary cores, the building blocks of planet formation, condense from the 
interstellar dust grains that are initially suspended in the gaseous protostellar disc. These first aggregate to form planetesimals 
of mass 10^^-10^'^ gm (e.g. Weidenschilling & Cuzzi 1993). Seed cores then grow up to an isolation mass, passing through a 
process of runnaway accretion (Lissauer & Stewart 1993). The duration of this growth phase is expected to be on the order 
of 10^ yr at 1 au (eg. see the reviews by Lin, Bryden, & Ida, 1999, and Papaloizou, Terquem & Nelson, 1999 and references 
therein). After isolation the orbits of cores do not overlap so that the system evolution slows significantly. Adopting a heavy 
element surface density of 6gmcm^^, characteristic of the minimum mass Solar nebula, the protoplanetary cores at isolation 
are thought to have roughly 0.03Me at 1 au from the Sun (Lissauer & Stewart 1993, Lin et al. 1999). For a surface density 
ten times larger this increases to IM^. 

The cores subsequently interact under their mutual gravitational attraction which can lead to collisions and the consequent 
production of more massive cores. If a critical mass scale ~ ISM,^ can be reached (e.g. Mizuno 1980, Bodenhiemer & Pollack 
1986), gas accretion can be initiated. Without further growth, the core can build up a substantial gaseous envelope which 
can contribute much of the final mass of the object as in the case of Jupiter and presumably the giant planets in extrasolar 
systems. The time-scale for the gas accretion process is estimated to be lO^'-lO^ yr (e.g. Bodenheimer & Pollack 1986). A 
potential problem is that the time required to build up a core with the critical mass may become very long if the cores become 
isolated and do not migrate (Lissauer & Stewart 1993). 

However, in order to discuss the evolution of protoplanetary cores, their gravitational interaction with the gas disc must 
be taken into account. This interaction produces wave excitation and angular momentum exchange (Lin & Papaloizou 1979) 
which can lead to orbital migration (Goldreich & Tremaine 1980). Low mass protoplanets interact linearly with the disc and 
undergo type I inward migration (Ward 1997b). Massive protoplanets interact nonlinearly and undergo type II migration 
(Lin & Papaloizou 1986). In this paper we shall be concerned with type I migration. In addition to orbital migration, 
eccentricity damping is produced (Goldreich & Tremaine 1980, Artymowicz 1993, 1994). Normally, migration calculations are 
undertaken assuming the eccentricity is smaller than the ratio of disc thickness to radius. Since the Solar nebula is generally 
believed to have been very thin, the latter is a small number, making the small eccentricity assumption very restrictive. 
Gravitational scattering amongst the cores acts to pump up their eccentricity until damping effects limit it. The possibility 
that the eccentricity is larger than the ratio of disc thickness to radius should be taken into account. We present here a new 
computation of the protoplanet-disc interaction in which it is calculated by summing over all resonances required to ensure 
validity, for eccentricities up to five times the ratio of disc thickness to radius, although in principle any value smaller than 
unity could be considered. 

We use our results to obtain an estimate for the equilibrium values of the mean eccentricity and vertical thickness for a 
system of protoplanets. This is obtained by balancing pumping through gravitational scattering with damping through tidal 
interaction with the nebula. The efi'ects of gravitational scattering and the details of the equilibrium are obtained by use of the 
Boltzmann Ti. theorem applied to the Fokker-Planck equation. The limitation of eccentricity because of the protoplanet-disc 
interaction is manifest in the restricted vertical extent of the protoplanet swarm. That is found to be always well confined 
within the gaseous disc, which in turn stabilises the collision rate and promotes rapid core agglomeration in comparison to 
the gas-free case. We investigate the consequences of this through numerical simulations. 

In Section 2 we calculate the dynamical relaxation time-scale for a system of equal mass protoplanetary cores by use of 
the Boltzmann Ti theorem applied to the Fokker Planck equation. In Section 3 we derive the eccentricity damping and orbital 
migration time-scales in the linear regime for a protoplanetary core that is fully embedded in a thin gaseous disc with surface 
density oc r~^'^ . In Section 4 we apply these results to find the vertical extent of a protoplanet swarm at equilibrium as a 
function of the nebula properties. In Section 5 we present numerical tests of our analysis using a direct summation N-body 
code in combination with an implementation of our nebula torque model. Specifically, we consider an ensemble of one hundred 
O.IM0 cores distributed interior to 1 au. 

The general finding is that, for characteristic protoplanetary disc parameters, because the eccentricity damping time-scale 
is significantly shorter than the migration time-scale, a quasi-equilibrium distribution is obtained in a lO'^-lO* yr time-scale 
at 1 au and which is otherwise proportional to the local disc radius. Inward migration occurs on the much longer migration 
time-scale and aids the accumulation of up to earth-mass cores in the inner regions of the disc on a 10^-10^ yr time-scale 
at ~ 1 au. However, these objects undergo rapid orbital migration towards the central star on the same time-scale, thus in 
this particular model, gas accretion onto a sufficiently massive core to make a giant planet within the disc lifetime requires a 
mechanism to halt the migration. One possibility is termination of the disc due to a magnetospheric cavity (Lin, Bodenhiemer 
& Richardson 1996). On the other hand if cores of the same initial mass are formed at much larger radii, they may survive 
for the disc lifetime without closely approaching the central star. In Section 6 we summarise and discuss these findings. For 
a first reading of the manuscript we suggest that the general reader moves directly to Section 4. 



2 ANALYSIS OF GRAVITATIONAL SCATTERING 

A system of many bodies interacting under their mutual gravitation can be described by the Fokker-Planck equation which 
we write in the form: 
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— Tcollifa) +Tgas{fa)- (1) 

Here fa denotes the phase space number density of bodies with mass m^. The operator giving evolution due to gravitational 
scattering is Fcoii, and that giving evolution due to interaction with the gaseous disc is Tgas- The latter combines the effects 
of orbital migration, and eccentricity and inclination damping. The derivative operator is taken following a particle orbit 
considering gravitational forces due to the central mass such that 

D d d ^ ^ d 

— = hv- V$ ■ — . (2) 

Where we refer to the central potential <1> — —GAI,/r, with the particle position and velocity vectors being denoted by r and 
V respectively. The central mass is M« and r = |r|. 

2.1 Local shearing sheet approximation 

As a result of the large relative velocities occuring between objects widely separated in radius, effective gravitational scatterings 
will occur on a radial length scale comparable to the vertical thickness of the planetesimal disc. The smallness of the latter 
quantity suggests use of the Goldreich & Lynden-Bell (1964) shearing sheet approximation. In this approximation we consider 
a uniformly rotating local Cartesian coordinate system in which the origin, located at some point of interest in circular orbit 
with radius/semimajor axis ro, corotates with the Keplerian angular velocity fi. The a::-axis points radially outwards and the 
y-axis points in the azimuthal direction while the z-axis points in the vertical direction. A linear expansion for the central 
potential is used such that 

V<1> = (nVo - 2r2^a;,0,n^2). (3) 

In this scheme we have for an axisymmetric disc with accordingly no dependence on y, 
D/a _ d,U , dfa , dfa on , /in , o ^o^/° o2 /a\ 

Dr = ^ + "^a^ + "^ar - 2^^"^ + ^^^^ + ^^^^^^^ ~ " 'a^' 

where v — {vx,Vy,Vz) and r — {x,y,z). The Liouville equation D/^/Dt — has from Jeans' theorem general solutions with 
fa being an arbitrary function of the integrals of the motion. Later we shall adopt such a solution, corresponding to an 
anisotropic Gaussian: 

fa-Ca^^Vy 2a|"2a2" 2a| j' 

Here the velocity dispersions ((t^, (Jy, (T2) are constant but such that Oy = (Jx/2. There is no constraint on a^. The velocity 
Uy = Vy -\- 3r2a;/2 is measured relative to the local circular velocity, Vy = —iQ.x/2, and we shall adopt local relative velocity 
vectors v = {vx,Uy,Vz)- The constant Ca is related to the spatial number density in the midplane, Ua, through 

Ca ~ I ,0/2 • (6) 

2.2 System evolution through gravitational scattering 

Encounters between planetesimals that occur without direct physical impacts tend to convert kinetic energy from shear into 
that of random motions, much as viscosity converts energy from shear into thermal motions in a gaseous disc. In the situation 
considered here evolution due to scattering occurs on a time-scale much longer than orbital. To investigate this phenomenon 
we use the form of Tcoii given by Binney and Tremaine (1987). This applies to a homogeneous system and so neglects rotation 
about the central mass. This should be reasonable so long as the time-scale associated with an encounter is short compared 
to . This in turn requires that a^jQ, > ru, where rn is the Hill radius f'.i^-^j^')^^'^ ■, appropriate to the characteristic mass 

Following Binney and Tremaine (1987), we write using the summation convention for repeated indices 

r..(/.) = -|-(^./.) + i|-(A.g). (7) 

Here 

A, =47rGMn(A)mi^ (8) 
and 

A, =47rG'ln(A)mi-^, (9) 

OViOVj 

with 
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(5, h) = l^j /.(v')|v - v'l d^v', j dV^ . (10) 

For A we take A — 3a^Ha / {'iGrria) , giving the ratio of maximum to minimum impact parameters as disc semi-thickness 
to impact parameter for a typical deflection (e.g. Binney & Tremaine 1987). The above formalism applies to one species of 
planetesimal with mass ma- Generalization to include a system of interacting planetesimals with different masses is straight- 
forward. However, for simplicity we shall assume all the planetesimals have equal mass in the analysis presented here. 



2.3 Growth of the velocity dispersion 

The effect of gravitational scattering is to cause the velocity dispersion to increase. This is most easily seen by formulating 
the Boltzmann Ti. theorem. This states that for a single mass species 

n = - j U\n{U)A\ (11) 

increases monotonically with time. Ti, which can be related to the entropy, can remain constant only for an isotropic Gaussian 
which cannot be attained here because of the form of particle orbits in the central potential (see equation ^). Using the 
distribution function given by (^) we obtain at the midplane (z = 0), assuming that is constant and for fixed ratio of 
velocity dispersion components, 

AH _ Sn^dCT^ 

At ^ At ' ^ ' 

Thus the velocity dispersion increases with time when ATL/At > 0. Here we make the assumption is constant and specialize 
to the midplane where 2 = 0. We comment that almost identical results are obtained if instead a vertical integration is done 
and the surface density is assumed constant. 

Evaluating AH /At, neglecting for the time being the effects of migration and damping, we obtain using (Q): 

— = -4.rGm.ln(A) J |,_,,| d v d v + ln(A) J ^ j-^^ j (^^-^j d vd v . (13) 

Note that the two terms on the right hand side of (jl^) give contributions of opposite sign. However, a net positive result 
is guaranteed by the Ti theorem which is why we follow that approach. The integrals in (113) are most easily performed by 
use of Fourier transforms and the convolution theorem with the result that at z = 0. 



^ = {4nfG^ml ln(A)(C.a.a,a.)2 J exp(-fc?a?) - 2j d^k. (14) 

When ay — az, the integral may be evaluated analytically. Using ( p^ we then obtain for the relaxation time, tR — ax/{Aax/At), 
1 8{ny^^G^mlln{A)n. 



tR 



(15) 



Here we have also specialized to the Keplerian case for which ay = ax/2. 

In general the vertical velocity dispersion would be affected by inclination damping. In the case of no damping we expect 
collisional effects to cause evolution towards isotropy such that — (ct^ -|- ay)/2, or a^/a^ = 0.79. Damping reduces this 
value, but as long as the rate of inclination damping is comparable to or less than that for eccentricity damping (as is expected 
to be the case from the results of Ward & Hahn 1994), this effect is not very significant (see below). To obtain a consistent 
estimate from ( |l5| ) we adopt the value of az/a^ ~ 0.5. We shall express the relaxation time in terms of the semi-thickness of 
the planetesimal disc Ha = az/^ and the planetesimal surface density Eq, — \/2tv H am aU a- Under these approximations we 
now have: 

(16) 

Thus our analysis is expected to break down when the disc thickness is smaller than the protoplanet Hill radius. Here we 
obtain: 

i.„.„3,„,A)^^(^)'a 

where Md ~ TrEarg gives an estimate of the protoplanet-disc mass within tq. Adopting values consistent with our numerical 
work: M, = IM©, ma = O.IM®, Md = lOM®, and Ha/ro = 0.01, we find for ro = 1 au that tR ~ 3000 yr. We comment that 
with ax/az = 2, Ha/ro — 0.01 corresponds to a mean eccentricity of 2y/2Ha/ro ~ 0.03. 

From this discussion we expect an equilibrium velocity distribution characterized as above, to have been obtained in a 
time ~ 10^ yr within 1 au, which is consistent with the numerical results presented below. The actual equilibrium levels of the 
velocity dispersion are obtained by balancing pumping due to gravitational scattering against damping due to tidal interaction 
with the nebula. We calculate the latter in the next section. 
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3 TIDAL TORQUES RESULTING FROM NEBULA INTERACTION 

Here we calculate the tidal interaction of a protoplanetary core with the gaseous disc. We assume the protoplanet mass is 
small enough that the disc response is linear. Thus we deal with migration of type I (Ward 1997b). 

We employ a simple two dimensional model disc for which the gas surface density a oc r"^''^, and the sound speed 
c oc r~^^^. The putative aspect ratio H/r is thus independent of radius and corotation torques are absent (e.g. Korycansky 
and Pollack 1993). The gas angular velocity Cl is then given by 

and as for Kepler's law this is proportional to r~^. 

The protoplanet or planetesimal of mass exerts a tidal potential 

, =, (19) 

y/r^ + 2rRcos{ip - (fa) + 

where {R,ipa) are the cylindrical coordinates of the protoplanet on its general Keplerian orbit. To prevent divergences we 
incorporate a softening parameter b. This can be thought of as taking into account the disc vertical thickness. 

The protoplanet potential is written as a double Fourier series (e.g. Goldreich and Tremaine 1978) in the form 

(oo oo \ 

>I'„,„expi[(n-m)a;t + mi^] J . (20) 
m— On— — oo / 

Here 5R denotes that the real part is to be taken and ui is the orbital frequency of the protoplanet taken to have orbital 
eccentricity e and semimajor axis a = ro. The Fourier coefficients are given by 

2-it/lj /■2tt 



*„„ = ——-—- r- / / '^cos[m.{ip-ipa)]difeyi^[-i{n-Tn)ut-irn^a\dt. (21) 

27r^(l + (5™,o) Jo Jo 

Where 5 is the Kronecker delta. Note that we do not include the indirect term in the above as the contribution due to this is 
much smaller than that due to terms with large m. When e = 0, '^nm ~ for all n 7^ 0. 

In the linear regime the torques due to each Fourier component may be evaluated separately and the results summed. 
The Fourier component (n, m) produces a tidal disturbance rotating with a pattern speed flp — (m — n)uj/m. In our model 
disc only Lindblad resonances are important and if they exist these occur at the two locations in the disc where 

2 2 

m^{Q,p — f2)^ — [(m — n)u) — mS7]^ — -\ — . (22) 

For m > 0, the outer Lindblad resonance with Q, <Q,p occurs where 



m(r2p - J7) = [(m - n)cj - mfi] = -I / H — . (23) 



The inner Lindblad resonance with Q, > Q,p occurs where 



m,(Q.p — n) = [(m ~ n)uj — mVl] = k? ~\ — . (24) 

Here the epicyclic frequency k = 0.. Note the term (P'rr? jr^ ^ normally neglected for small m, leads to an accumulation of 
resonances at the two locations where (tJ — fl) = ±c/r as m — > 00 for finite n. These locations are at a distance 2iif/3 inside 
and outside the semimajor axis of the orbit. There are no closer resonances for large m. This leads to the well known torque 
cut off for m ^ r jH (Goldreich & Tremaine 1980, Artymowicz 1993) which results in convergence of torque sums in the 
case of a circular protoplanet orbit (n = 0), even when gravitational softening is omitted. However, for general n, resonances 
may approach the coorbital radius and so results become more sensitive to the softening parameter. For our disc model the 
resonant locations occur where 

^ ("^ - (25) 



Thus a strictly coorbital resonance occurs if n = ±-^1 + {H'^m? /r^). Even when |n| = 1 two resonances occur at locations 
satisfying {Q. — uj) ^ {H/r)uj when m ~ r / H. These have been found to lead to rapid eccentricity damping (e.g. Artymowicz 
1994). Up to now terms with jn| > 1 have not been considered, which is reasonable for very small e. Here we consider all n 
necessary for convergence. Once e > H/r, significantly larger values of n are required and these may produce more closely 
coorbital resonant effects. 

We comment that for n = 0, the inner Lindblad resonances are all located interior to the semimajor axis while the outer 
Lindblad resonances are all located exterior to it. This results in the inner disc unambiguously providing a positive torque and 
the outer disc unambiguously providing a negative torque. However, when n > 0, inner/outer Lindblad resonances may fall 
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Figure 1. The magnitude of the cumulative torque in arbitrary units as a function of m. The net torque contributions arising from 
inner Lindblad resonances are plotted with solid lines and those arising from outer Linblad resonances are plotted with dashed lines. In 
the left-hand panel we show the magnitudes of the net torque contributions for very small e and H/r = 0.07. In this case only n = 
contributes for each m. The outer Linblad resonances predominate leading to net inward migration. In the right-hand panel we show the 
magnitudes of the net torque contributions for e = 2H/r and H/r = 0.07. In this case the torque at each m has all of the significant 
contributions from each value of n included. In this case, unlike when e is small, the inner Lindblad resonances predominate, leading to 
net outward migration. Note that these resonances are not all located in the inner disc. 



external/internal to the orbital semimajor axis. This means that the interior or exterior discs may produce torques of either 
sign on an orbit with significant eccentricity. Physically, a protoplanet at apocentre may rotate more slowly than coorbital 
disc material. If e is large enough compared to H/r the exterior disc can then exert a positive torque. We apparently find 
that this phenomenon may lead to a net torque reversal for sufficiently eccentric orbits. 



3.1 Rate of orbital evolution 



To evaluate the rate of change of protoplanet angular momentum arising from the (n, m) Fourier component, we use the 
Goldreich and Tremaine (1978) torque formula as modified by Artymowicz (1993) (see also Ward 1997b) in the form 



dJ„,; 



At 



dr 



+ 



2m^(r2- r2p)*„,, 



H-4m2cV(r2!:!2 



(26) 



Here S — 1 for an inner Lindblad resonance that transfers angular momentum to the planet, S = —1 for an outer Lindblad 
resonance that removes it from the planet and the right hand side of (^) is to be evaluated at the resonant location. The 
associated rate of change of orbital energy is then related to the rate of change of angular momentum through: 

dE 

,m 



dt 



(27) 



where we have added the subscripts [n, m) to the pattern speed Q.p to denote that it is associated with the corresponding 
Fourier component. 

By summing over (n, m) we find the total rate of change of angular momentum dJ/dt and we define the orbital migration 
time tm such that 



tn 



J 



(dJ/dt) ' 



(28) 



where J is the angular momentum of the protoplanet. Migration here is thus defined in terms of the total torque exerted on 
the orbit. The migration time so defined is positive when the total torque is negative. 

Similarly we find the total rate of change of orbital energy and then the rate of change of eccentricity using 



dE 



J \dt 



( 



dJ 



1 dJ 



J dt 

The eccentricity evolution time-scale is then 
e 



de 



(l_e2)3/2 



te 



\de/dt\ ■ 



(29) 



(30) 



When, as is always found to be the case here, de/dt < 0, t^ gives the circularisation time-scale. 

We have performed calculations of tm and te taking into account as many resonances as required for convergence. We 
comment that we did not make any large m asymptotic approximation as this was found to affect the results significantly. 
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We illustrate the contributions to the torque for very small e for which only n — and n = 1 need to be considered in the 
left-hand panel of Figure 1. The calculation illustrated is for H/r = 0.07. In this case the contribution from the outer Lindblad 
resonances that are located in the outer disc dominates and the migration is inwards. The eccentricity is always damped. 

We found that departures from the small e case occured once e > H/r. The sign of the torque actually reverses once e 
exceeds ~ l.lH/r. Significant changes occur because for such values the radial excursion causes the protoplanet to cross the 
circular orbit resonances a distance 2H/3 away from the orbital semimajor axis. 

In the right-hand panel of Figure 1 we show the torque contributions from the inner and outer Lindblad resonances for 
H/r = 0.07 and e — 2H/r. In this case the inner Lindblad resonances dominate giving a net positive torque. But these are 
not always located at radii interior to the orbital semimajor axis. The eccentricity is always damped, but note that te, oc 
for e > H/r, so that the more eccentric orbits feel very much weaker damping than the near circular ones. 

We found that the dependence on the gas-disc thickness was such that tm oc (H/ro)'^, and te oc (i//ro)'*, with H now 
being evaluated at r = tq. This scaling, together with the coorbital resonance contributions to t^, result in being much 
shorter than tm. Accordingly we expect an equilibrium eccentricity distribution to be set up in a system of gravitationally 
interacting protoplanets embedded in a gas disc, assuming that they are massive enough so that tidal torques dominate over 
those due to gas drag. After this equilibrium is set up, longer time-scale evolution is expected to occur due to orbital migration 
and physical collisions. 

We also note that the torque results were dependent on the softening parameter used with tm and te, scaling as b^'"^^ and 
b^'^ respectively, for b in the range 0.4H-H. This leads to an uncertainty reflecting the inadequacy of using a two dimensional 
flat disc model as the softening would occur naturally in a fully three dimensional treatment. Fortunately our simulations of 
protoplanet discs are insensitive to this issue. Also as we find the eccentricities are for the most part limited to small values, 
making the results for this case the most important for the work presented here. For small e our results are broadly consistent 
with those of Ward (1997b) and Artymowicz (1993, 1994). From our calculations we found the approximate fits to tm and te 
given by 




Here the gas disc has a mass AIgd contained within 5au. At other radii the mass of the gas disc is assumed to scale as Tq , 
as implied by the model surface density. The Jovian mass is Mj. We note that the factor due to softening, fs = (^^), is 
hereinafter taken as unity throughout, coressponding to 6 = O.iH. For H = 0.07ro, Mod = 20Mj, and rria = O.lAf^, we 
obtain te ~ 10'^ yr and tm ~ 10^ yr, in the limit of e ^ H/r at 1 au. 



3.2 The contribution from inclination damping 

In estimating the damping forces we return to the Boltzmann Tl theorem and consider the contribution to dTl/dt one obtains 
from r gas if a) in equation (|l]). We suppose this can be represented through the action of a body force F per unit mass such 
that 

rgas{.U) = ~-^{FU). (33) 

Here if F = —{2vx/te,0,2vz/ti), the eccentricity and inclination damping times will be te and ti, respectively. As migration 
occurs on a longer time-scale, we shall neglect it when considering the velocity dispersion equilibrium. The damping force F 
gives a contribution to the rate of increase of Tt in the midplane: 

Thus, as long as ti ^ te, the equilibrium condition dTC/dt — can be approximated by equating te to the relaxation time tu, 
see equation (|l7|). 



4 VELOCITY DISPERSION EQUILIBRIUM 

As we have seen in Section ^ the effect of gravitational scattering not involving direct collisions is to increase velocity dispersion 
on a time-scale t^. An equilibrium is attained by balancing this rate of increase against damping due to tidal interaction of 
the protoplanet with the gaseous disc. As described in the previous Section we can approximate damping effects by ignoring 
inclination damping, occurring on the time-scale ti, and only consider eccentricity damping, occurring on the time-scale te, 
provided tei^ti. 

The time-scale te is obtained from the protoplanet-disc interactions calculated above. As this time-scale is significantly 
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shorter than the migration time-scale tm, we expect that a quasi-equihbrium is set up before significant migration occurs and 
that this then drives the evolution of the system, which will occur under conditions of quasi-equilibrium. To calculate the 
quasi-equilibrium we use our fit to te given by ( p2|) above and recall our result on the dynamical relaxation time-scale for a 
central Solar mass point potential: 



3/2 

tR 



In (A) Momo 



where A is given by equation (jl^. Equating te from equation and ta, assuming e < H/ro gives the equilibrium semi- 
thickness of the protoplanet distribution Ha in terms of the gas disc semi-thickness H: 



/ rn \ -1/8 , ^ 

(36) 



^=0.6[ln(A)]^/^f#^)" ff^ 
H L V yj \Mgd) Vlau 

Ha/H depends only very weakly on A and therefore (|6|) is essentially independent of the protoplanet mass m^, and depends 
only weakly on other parameters. For rria = O.IMq, protoplanet disc mass Md = lOM^, Ha/ro = 0.01, and gas disc mass 
within 5au Mod = 2QMj, we determine Ha = 0.15i/ at 1 au. Hence the protoplanet swarm is expected to remain thin 
and confined within the gaseous nebula. The characteristic eccentricity of the equilibrium (see Section 2) is then roughly 
2y/2Ha/r ~ QAH/r. 



5 NUMERICAL MODELLING 

In order to investigate the consequences of our approximate analysis, we have developed a three dimensional direct summation 
N-body code to model the gravitational dynamics of a system of protoplanetary cores. For the time integration we employ a 
fifth order accurate Runge-Kutta-Fehlberg scheme with time step control through the estimated error in velocity magnitude 
(Press et al. 1992). The scheme is not strictly conservative with the result being orbital decay in the circular orbit two body 
problem. The accuracy tolerance - used in calculating the computational timestep - was calibrated through experiments with 
a one Jupiter mass body in circular orbit about the Sun at 0.1 au. We set the accuracy tolerance such that the orbiting body 
suffered 0.01% orbital decay in 10^ yr, being a much weaker effect than orbital migration due to nebula interaction. 

To incorporate the effects of eccentricity and inclination damping and orbital migration owing to tidal interaction with 
a gaseous disc we implemented the following expressions as contributions to the total acceleration for each particle: 

V 

^damp = 2-^ 2 . ^ 7 — ^ — ' {^^) 

T ''i 

Here k is the unit vector in the vertical direction. As long as ti is not significantly less than t^, tests confirm that inclination 
damping does not change the velocity dispersion equilibrium significantly (see below). Hence the root-mean-square velocity 
dispersion of the particle distribution is essentially controlled by the quasi-equilibrium eccentricity attained by balancing 
pumping through gravitational scattering with damping owing to tidal interaction with the gaseous disc. 

In all our numerical tests we consider 100 equal mass particles of O.IM® each, and a centralised Solar mass. The particle 
distribution is initialised in the gas-disc midplane such that the surface density is proportional to r~'^/^. The positions are 
otherwise chosen randomly in the radial range 0.3-1 au such that the separations are greater than twenty Hill radii. The initial 
velocity components in the midplane correspond to Keplerian circular motion about the central object and vertical velocities 
are defined as iv^, where i is chosen at random for each particle to lie between ±0.01, and is the Keplerian circular velocity. 

In a real protostellar accretion disc we expect the mass of the disc at a fixed radius to decrease with time on the global 
viscous time-scale, which at 1 au is typically ~ 10'* yr (Papaloizou & Terquem 1999). We note that te oc H'^/Mgd and 
tm oc H^/Mgd, which tend to increase with diminishing gas-disc mass. However, this increase is compensated for by the 
decreasing H found in nebula evolution models (Papaloizou & Terquem 1999). Here we consider a fixed value Mod ~ 20Mj 
which corresponds to an early stage in the life of the nebula at an age ~ 10*1-10^ yr, in our numerical models. The gas disc 
then has a mass ~ O.IMq within a radius 100 au. However, our calculations of the equilibrium eccentricity distribution (see 
below) can be scaled to a lower gas-disc mass with smaller vertical scale height by scaling H cx Mq^. They can also be 
scaled to apply to larger radii by the scaling transformation: (r Xr,t — > X^'^t, Mod MgdX~^^^). Thus all radii can be 
multiplied by 9 provided all times are multiplied by 27 and the disc mass within 5 au is reduced by a factor of 3. 

In some of our numerical models we allow the masses to collide and agglomerate. We simulate this process by replacing 
closely approaching particle pairs with a single mass occupying the centre of mass position of the original pair: the new masses 
are set up to conserve the mass and linear momentum of the original particle pairs. Binary agglomeration occurs when the 
separation is less than the sum of the radii associated with each component. For this purpose each mass is assigned a radius 
that gives a mean density equal to the lunar value. This approach assumes that the relative motion of the colliding pair is 
completely dissipated and that the liberated kinetic energy does not exceed the total binding energy of the putative bodies. 
For the relatively large mass scales considered here the latter is expected to be realistic (Lecar & Aarseth 1986) . Furthermore, 
modelling of coUisional fragmentation shows that > 90% of the total mass remains in the merged core and the existence of 
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time (yr) time (yr) 

Figure 2. Time evolution of the mean eccentricity (uppermost curves) and the mean value of \z\/r (lowermost curves). The left-hand 
panel is for a model without nebula torques included in the force calculation. The right-hand panel shows data for a run with eccentricity 
damping only (solid lines) and a run with both eccentricity damping and inclination damping (dotted lines). 




a (au) a (au) 

Figure 3. The eccentricity and vertical distributions of particles in the model having H/r = 0.1 and no migration or inclination damping, 
In the right-hand panel the dashed line refers to the scale height of the putative gas disc. 



remnant ejecta is not significant in the subsequent evolution of the planetary system (Alexander & Agnor 1998), and thus 
can be ignored. 

We note that if agglomeration is included in the simulations as described above, and the A transformation is applied, 
the physical radius of cores is also scaled. This has the affect that for example if A > 1 the agglomeration time-scale is 
underestimated. None the less the transformation still serves to indicate how evolutionary time-scales can be increased if 
initial cores of the same mass are moved to larger radii with a correspondingly reduced disc mass. 



5.1 Equilibrium eccentricity 

In Figure 2 we compare the time evolution of the mean eccentricity and the mean value of \z\/r, which is used to estimate 
the characteristic value of Ha/r. We do this for models with and without damping, but for which we did not include orbital 
migration or agglomeration. We recall that the dynamical relaxation time-scale estimated above is ta ~ 1400 yr, at the mean 
radius. For the model that includes gas with H/r = 0.1 we find ~ 6800 yr, at the mean radius for e <C H/r from equation 

(§)• 

Without nebula interaction the model cores show eccentricity and vertical thickness increasing monotonically with time. 
The model that only uses eccentricity damping shows that a quasi-equilibrium in the mean eccentricity is attained after 
~ 4000 yr. The average eccentricity is ~ 0.3H/r and we infer Ha/r ~ 0.3H/r, broadly in agreement with our analytical 
estimates outlined in the previous Section. 

Comparison of the eccentricity-damped model with a similar model that also includes inclination damping with ti = te 
shows that the time-scale and level of the equilibrium eccentricity is essentially the same as before. In the former model, the 
equilibrium thickness is established after a time ~ 4000 yr. In the latter model, the equilibrium value is established after a 
time ~ 2000 yr, at a level of about 40% less than for the former model. This value is minimal given the expectation of ti >te. 
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Figure 4. Time evolution of the mean eccentricity and the mean value of |2|/r, for models with damping, migration and agglomeration. 
The left-hand panel shows data for the H/r = 0.07 runs with eccentricity damping only (solid lines) and with inclination damping 
included (dotted lines). In both cases the mean eccentricity has the higher value. The right-hand panel shows data for the H/r = 0.05 
model (uppermost curves, for which the higher value corresponds to the mean eccentricity), and also the H/r = 0.03 model (lowermost 
curves, for which the mean eccentricity initially has the higher value). 

These calculations assumed Mgd = 20Mj, but can be compared to other gas-disc masses by scaling H oc Mq^; for example 
the minimum mass Solar nebula with Mgd ~ 2Mj has the same equilibrium as the model used here when H/ro — 0.06. Thus 
the aspect ratio of the protoplanet swarm relative to the gas-disc aspect ratio increases as the gas-disc mass is reduced. 

In Figure 3 we give the eccentricity and vertical position distributions at a time of 10* yr, for the model without inclination 
damping. The maximum eccentricity is ~ 2/3H/r. The low level of eccentricity and inclination in all cases with damping puts 
the protoplanet distribution well inside the putative disc gas environment, as demonstrated in the right-hand panel of Figure 
3. This situation promotes physical collisions in a system of bodies with finite size, which we consider next. 

5.2 Core agglomeration models 

We present models with H/r = 0.07, 0.05 and 0.03, employing both eccentricity damping and orbital migration. Again, 
Mgd ~ 20A/j. For comparison we repeated the model with H/r — 0.07, but augmented with inclination damping such that 
ti — te- Additionally, there is one model without damping or migration. For computational convenience the radius of the 
central object is taken as lOi?0 in all cases with orbital migration. First we consider the previously discussed equilibrium in 
agglomerating simulations, and then protoplanetary growth and orbital migration. 

5. 2. 1 Equilibrium 

Figure 4 shows the average eccentricity and the mean value of \z\/r for all our models having migration and damping. The 
runs with H/r = 0.07 are qualitatively similar to the analogous runs with H/r = 0.1 of the previous Section, with a value 
of ~ 0.3H/r for the average eccentricity and ~ 0.2H/r for the semi-thickness to radius ratio, in the case without inclination 
damping. Again, inclination damping has little affect on the eccentricity equilibrium, and a moderate affect on the particle-disc 
thickness, being reduced by ~ 40% in this case. The model having intermediate gas-disc thickness H/r = 0.05 gives ~ 0.2H/r 
for the average eccentricity, and for Ha/r, but we note that the averages show a long time-scale decline. This is due to the 
diminishing particle number, which is a stronger effect in the thinner gas-disc models than for the thicker gas-disc model (see 
below). We note also that in this model the relative mean eccentricity has dropped in comparison the previous cases, and it 
is about a half of the predicted value. However, incorporation of migration and agglomeration appears not to strongly affect 
the establishment of the predicted velocity dispersion equilibrium in these cases. 

The thin disc model with H/r = 0.03 shows an equilibrium semi-thickness to radius ratio ~ 0.2H/r, but the average 
eccentricity is much smaller than for the thicker disc models at ~ O.lH/r, although this has a value ~ 0.3H/r for the first 
2000 yr. In this model the particle-disc thickness is comparable to the mean Hill radius. Consequently the curvature of the 
Keplerian orbits makes scatterings less efficient at pumping eccentricity than predicted and so the average eccentricity seeks a 
lower equilibrium value. Hence it is likely that the intermediate gas-disc thickness model is close to being marginally applicable 
to our analysis, and the agreement on the particle-disc thickness found in the thin gas-disc model is coincidental since the 
initial mean Hill radius is ~ 0.1-ff/r. In our models, in order to compare behaviour at small H/r < 0.05 with the analytical 
predictions, we would require less massive planetesimals. However, since our analytical results are relatively insensitive to m^, 
similar qualitative behaviour would be expected. 
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Figure 5. The total particle number in the simulations versus time. The data for the model without nebula interaction is plotted with 
a dashed line, the data for the model with H/r = 0.07 and inclination damping included is plotted with a dotted line. From highest to 
lowest final value the solid lines are for models with eccentricity damping and orbital migration having H/r = 0.07, 0.05, 0.03. 

5.2.2 Protoplanetary growth and orbital migration 

In Figure 5 we show the total particle number versus time for models with and without nebula interaction. The former 
calculations include both damping and orbital migration. The agglomeration time-scale for cores is reduced by including 
nebula interaction, owing to vertical confinement of the particle distribution, limited radial velocity dispersion and orbit 
crossing due to migration. The rate of agglomeration for the model with ti — t^ and H/r = 0.07 is of similar order to that 
with H/r — 0.05 and t^ only, being only about 20% faster than that with H/r = 0.07 and te only. Hence inclination damping 
does not very significantly enhance the agglomeration time-scale provided ti ^ te. In all cases with gas-disc models included 
the accumulation rate is much faster than for the gas-free model. We note that the model with H/r — 0.07 has a migration 
rate that is comparable to that expected for the minimum mass Solar nebula with Mod ~ 2Mj but having H/r — 0.02. 

In Figure 6 we compare the particle trajectories for the models with eccentricity damping and migration only and the case 
without nebula interaction. In all the cases with nebula interaction included we find that several ~ O.S-IM® cores build up 
and undergo increasingly fast migration to the central regions such that the time-scale to build up these masses is comarable 
to their net migration time-scale. The system thus appears to tend to the situation in which the remaining low mass cores 
are sufficiently widely spaced that scatterings are no longer effective in raising the eccentricities against damping. However, 
as long as the core masses are not uniform, differential migration of cores allows further accumulation and rapid migration 
until uniformity is achieved, or all matter is delivered to the central star. If the former situation occurs, since tm increases 
linearly with radius, there can be no further mass accumulation, as the cores' orbital separations do not decrease with time. 



6 DISCUSSION 

We have derived eccentricity damping and orbital migration time-scales for protoplanetary cores fully embedded in a model 
gaseous nebula disc with E oc r"'^'^. This has been done for e significantly larger than H/r, taking into account for the first 
time all effective contributing resonances (n,m). We find that for 1 > e 3> H/r circular orbit calculations are modified in a 
way indicated by the replacements: te [e/{H/r)]^te and tm —e/{H/r)tm. 

This behaviour is in accord with the simple impulse model for tidal interaction (Lin & Papaloizou 1979). According to 
this the coupling between the orbit and disc would be expected to be weaker for eccentric orbits owing to the larger relative 
velocity of the protoplanet and local disc material. The weaker gravitational interaction explains the generally longer time- 
scales obtained in comparison to the circular orbit case. We also find reversal of the direction of migration, defined according 
to whether the orbital angular momentum decreases or increases, is possible for sufficiently eccentric orbits (e <: l.lH/r). 
This can be thought of as a consequence of the protoplanet spending most time interacting with the disc at apocentre, where 
the relative velocity of the disc matter gives a predominating contribution to net migration in the outward direction. 
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Figure 6. Particle semimajor axis versus time for the models without inclination damping, plus the model without nebula torques. Note 
that the data is not uniformly spaced in time. 



We note also that the circularization time-scale can be significantly increased when the orbit has significant eccentricity. 
Application of ( |32[ ) indicates that for Mqd ~ 2Mj and H/r — 0.07, a critical mass core of ISM® has a circularization time 
of ~ 1.6 X 10^ yr at 10 au. If e = 0.35 this is increased to 10^ yr, which is comparable to the local nebula accretion time- 
scale. Thus standard inward migration can be prevented if the eccentricity can be maintained against decay by some process 
such as gravitational interaction with either larger objects and/or or a distribution of planetesimals with enough total mass. 
Accordingly gas accretion onto a critical mass core in an orbit with modest eccentricity should be studied. In this context 
we note that the notion of a clean annular gap being opened in the disc by the action protoplanetary tides breaks down for 
eccentric orbits because resonances giving positive and negative torques exist at radii both less than and greater than the 
semimajor axis. In the small eccentricity case the positive and negative torques are divided separately between the opposing 
sides of the orbit (Ward 1997b). 

We calculated the dynamical relaxation time-scale for a system of cores and deduced the equilibrium eccentricity set up 
when there is a balance between pumping through gravitational scattering and damping through nebula interaction. We find 
that the vertical distribution of the protoplanet swarm generally remains well confined within the gaseous envelope of the 
nebula. This supports our use of a two dimensional analytical model for the tidal interaction. 

We determined analytically, and confirmed empirically, that the mean eccentricity and the thickness to radius ratio for 
the protoplanet swarm are roughly 20-30% of the gas-disc aspect ratio, provided the characteristic size of a protoplanet 's 
Hill sphere is smaller than the latter. But we note that there is a weak dependence on gas-disc and particle-disc mass, 
and also radius. If the characteristic Hill radius is larger than or comparable to the thickness to radius ratio then lower 
equilibrium values of thickness and eccentricity occur. The vertical confinement of the cores enhances the collision frequency 
and consequently promotes more rapid agglomeration into larger objects. 



6.1 Protoplanetary growth 

With the aid of a numerical code we showed that the presence of a gaseous disc can dramatically affect the coUisional evolution 
of a system of gravitating cores. We found that several cores increased in mass by a factor ~ 10 in only ~ 10* yr, and underwent 
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migration to the central star on the same time-scale, for a gas-disc ~ 10 times more massive than the minumum mass Solar 
nebula. Our scaling for tm is such that the migration time-scale for a protoplanet in a Mqd = 20Afj gas disc is the same 
as for the same protoplanet in the Mgd ~ 2M.j minimum mass Solar nebula gas disc, but with a moderately reduced H/r. 
Additionally, the migration rate of our adopted O.IMq initial core isolation mass at lau scales to the same value for a IM^ 
protoplanet at 10 au. Starting from an ensemble of equal mass cores embedded in the disc model with E oc r~^'^ , as described 
above, implies that rapid gas accretion onto a massive core cannot occur before it plunges into the star. Tanaka & Ida (1999) 
reach similar conclusions by considering the growth rate of a migrating IA/q protoplanet through a field of much smaller 
planetesimals. By extrapolation of the growth rate, they find that for the minimum mass Solar nebula, protoplanets could 
only build up to a few earth masses before being consumed by the stellar envelope. 

Thus the inward migration of at least one core must be halted short of the stellar envelope by some means at the initial 
stage of clearing, if a gas giant is to be formed at all. If a halting mechanism, such as the presence of a magnetospheric cavity, 
operates then other cores would be delivered to the first one owing to orbital migration. Providing the migration is not too 
rapid the inner body could feed on the inflowing objects (Ward 1997a, Papaloizou & Terquem 1999), presumably giving a 
critical core in only ~ 10" yr. But we note that no halting mechanism is known to apply to Jupiter at its present location of 
5 au. Although it may be possible that the cores of the outer planets were formed at much larger disc radii with correspondingly 
longer migration time-scales, and were possibly much closer together in the very early stages of their formation, which could 
also allow for the mutual excitation of their eccentricities and the consequent stalling of inward orbital migration. It could also 
be possible that planetesimal scatterings by the cores could help to maintain their eccentricities (Hahn & Malhotra 1999). 

Vertical confinement of the protoplanet swarm through inclination damping helps to reduce the agglomeration time-scale 
without increasing the migration rate. However, for inclination damping time-scales greater than the eccentricity damping 
time-scale this is not a very significant effect. Even when the zero thickness limit is approached, we still have the problem of 
type I orbital migration causing very rapid inward migration of the most massive cores. 
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